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Abstract In a recent paper the authors studied numerically the hydrogen ground 
state in stochastic electrodynamics (SED) within the the non-relativistic approx¬ 
imation. In quantum theory the leading non-relativistic corrections to the ground 
state energy dominate the Lamb shift related to the photon cloud that should cause 
the quantum-like behaviour of SED. The present work takes these corrections into 
account in the numerical modelling. It is found that they have little effect; the 
self-ionisation that occurs without them remains present. It is speculated that the 
point-charge approximation for the electron is the cause of the failure. 

Keywords Stochastic Electrodynamics • hydrogen ground state • relativistic 
corrections • simulations • OpenCL 

PACS 11.10 • 05.20 • 05.30 • 03.65 


1 Introduction 

Stochastic Electrodynamics (SED) has been put forward as a candidate for a 
hidden variables theory that underlies quantum mechanics, see the books [1,2]. 
Such theories are not ruled out by Bell’s theorem, since that has an irreparable 
fatality, the contextuality loophole [3]. 

The derivation of the Schrodinger equation from SED has been reported [4, 
5] including “truly quantum” properties like the photoelectric effect. A study of 
radiative corrections opens up the contact with quantum electrodynamics [6]. A 
parodox between initial classical and final quantum behaviour in scattering exper¬ 
iments is explained [7]. Finally, it is realised that the energy throughput of SED to 
keep matter stable (energy from stochastic fields in, raditation out) can be seen as 
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an arrow of time, called the subquantum arrow of time, which is more fundamental 
than the entropic and cosmological arrows of time [8]. 

SED is known to work well for harmonic (linear) problems, see e.g. [1]. A 
fundamental test for nonlinear problems is the ground state of the hydrogen atom. 
While Claverie and Soto put forward that the motion is non-recurrent, and can 
not have a stationary distribution [9], Puthoff puts forward a comparison between 
energy gain and loss terms, concluding that it should be stable [10]. 

The theory has been tested numerically on the hydrogen ground state in a 2-d 
approximation by Cole and Zou 2003, who observed an encouraging agreement 
with the result from the quantum theory [11]. New simulations were carried out 
recently by us (ref. [12], to be called NLl) without the 2-d approximation and 
including much more computational detail and power. However, it was observed 
that in all runs and all modellings of the system, there occurred self-ionisation. 

In quantum theory the contributions from the leading relativistic corrections 
go as 1/c^, hence as with a = 1/137 the hne-structure constant. The SED 
correction to the ground state energy is related to the Lamb shift, which is of 
order a^loga. The relativistic corrections, mechanical terms, are formally more 
important than the SED effects, whereas in SED theory they should only act 
as small contributions to the non-relativistic Hamiltonian. To see whether this 
paradox plays a role in the final stability question, their role has to be analysed 
numerically. For consistency, the stochastic electromagnetic vector potential has 
to be calculated to relative order as well, that is to say, it has to include spatial 
terms of order r and . 

We take up this challenge and compare the results with the conjecture for the 
phase space density of the ground state [13]. This theory will be recalled in section 
2. Next the electromagnetic fields, Gaussian sums over 3d fc-space, are replaced by 
Gaussian sums over frequency that reproduce the same correlation functions. This 
essential simplification allows us to carry out simulation of the dynamics, which 
we program in OpenCL. The setup is explained in section 4 and the results are 
reported in section 5. In section 6 we close with a discussion. 


2 Relativistic corrections in the hydrogen problem 

The spin-orbit interaction reads for an electron of mass m and charge —e in the 
field of a nucleus with charge Ze (e > 0), with Bohr magneton fis = eh/2m in SI 
units. 


Hso 


d^L s 

Timec^r dr 


The relativistic Hamiltonian including the Darwin term reads 


( 1 ) 


Rrel = 


Ze^ Ze^ L-S 
dTreor STreom'^c^ 


( 2 ) 


The Bohr units 
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ry ’ ^0 rjo 2 2 ’ \'^/ 

Zoivac ijjQ Z^a^mc^ 

allow to introduce dimensionless variables by r —^ agr and p — ^ (mao/To)p, 
implying that L = r x p —^ hL, and to take consistently S —)• ?iS. Keeping 
|S| = one arrives at 


mc^ 


yi + Z2a2p2 _ 


/72 2 
Z a 


, ZWL-S , ^4 , 

+ ^—+Z^^5{r) 


(4) 


So the dimensionless nonrelativistic Hamiltonian = {Hre\—'m(?)/a^m(? picks 
up the leading relativistic corrections 


u 1 2 1 1^72 24, 2 2 7r ^ ^ 

^=2P-r-8^^P +^—+ ^« 2^W- 

From r = dpH and p = —drH the dynamics reads for r ^ 0 


(5) 


. 2 s ^ a S X r 

r={l-^P )P+— - 

_ r , Z^a^ S X p , dZ^a^ L • S ^ 

P - 2 ~ 2 


while the spin progresses as 


(6) 

(7) 


Z^a^ L X S _ Z^a^ (pr- rp)-S 

2 y .3 2 


(8) 


It is a standard exercise to demonstrate that the Hamiltonian H and the total 
angular momentum J = L + S are conserved, here up to terms of order a^. 

An external electromagnetic field is added by the minimal substitution p —> 
p = p + /3A(Zar, t), with /3 given below and the spatial scale factor Za expressing 
the ratio of the Bohr radius to the wavelength of a photon with Bohr energy a^mc^. 
From Eq. (8) it is conhrmed that remains conserved, as desired. 

Stochastic Electrodynamics leads to a specihc stochastic electric and magnetic 
held, as well as to an r' damping term, see NLl and Ref. [1]. When we neglect 
terms of order a^^^ higher, the contribution Vip to (7) and the time derivative 
of (6) lead to the Abraham-Lorentz or Brafford-Marshall equation for a particle 
with spin. 


f = - ^ - ^(E + f X B) + /3^ r + Z^a^ 


r^r + 2f r f 


2r^ 


Z^a\^ . r-r„ S-L , 

+ X *•-3^8 xr + 3^r). 


(9) 


where both the huctuations and the damping involve the small parameter 


f Za^/^ = --, 

3 1964.71 


1 


dTreohc 137.036 


(10) 
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(In order to have /3 also as prefactor of E, we absorb a factor ^3/2 in A, B 
and E.) We employ the standard rule dA/dt = dA/dt + r- V A and definitions 
E = —dtA, B = V X A. Here E is the transversal part of the electric field, next 
to the longitudinal part, viz., the Coulomb force —r/r^. 


3 Representation of the stochastic fields 

In the Coulomb gauge the SED vector potential of a cube of volume E is a sum 
of plane waves with random coefficients. 


A(r,t) = ^4/-—rye [ AkA sin(k • r - ojkt) 

tt V 2eo^^kE 

+ SkA cos(k • r - ojkt)] (11) 

The wave vector components ka = 2'nnalV^^'^ involve integer Ua = —oo, • • • , oo, 
(o = 1, 2, 3). The £kA with A = 1, 2 are polarisation vectors. The AkA and SkA are 
independent random Gaussian variables with average zero and unit variance. For 
each (k. A) term the average energy d^r(^(E^) + 2 ^(B^)) = yliixik exp(—ojkTc) 
is equal to the photon zero point energy combined with an exponential cutoff; we 
choose Tc = ‘hjrat?, the Compton time of the electron. 

3.0.1 Noise correlators 

Because it is not often listed in text books, we start with deriving the autocorre¬ 
lator of the A held C^(r, t; q, s) = (A(r, t)A(q, s)). It is translationally invariant 
in space and time, viz. Czi(r, t; q, s) = Czi(r — q, t — s) with 


CA{r,t) = ^ ‘^‘‘^“cos(k-r-a;kI). (12) 

k 

Symmetry considerations yield the decomposition 


Ca = Cq! — Ciff, trCA = 3Co — Cl, trCzi • ff = Co — Ci. (13) 
Let us introduce Cs and Cp as 

Cs = itrC, Cp = itrff-C, (14) 

which determine 

Co = Cs- Cp Ci = Cs- 3Cp. (15) 


Integration over angles k yields 
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Cs{r,t) 
Cp{r,t) 
where a = c{t 


h 


h 


47r2eoC 

n r 


■5ft 


1 


47r2eoc 
— ir). The functions 


,, , _T-cfcSinfer , , 

dfc ke —;-cos ckt = — , „ , 

kr 

' , , --rh sin kr — kr cos kr , , 1 / ^ 

dfc fee - cos ckt =^{Fi- 


{kr 


(16) 

F2m 


Ti = 

n 

fOO 

/ dfc ke 

'to 

dTT^eoc . 

F2 = 

n 

fOO 

/ dfc ke 

'to 

dvr^eQC . 


— Tck 


— Tck 


1 — cos kr 
kr — sin kr 


k^r 


cos ckt, 
cos ckt, 


(18) 


can be found from Cs since drFi = d^irF^) = rCs{r), while _Fi ~ F 2 ~ r^ for 
r —)■ 0. They read 


Fi = 
F 2 = 


n 


Stt^coc 

n 


3ft log 


2 2 

a — r 


3ft 


2 2 

a — r 


a cr + r 

- log-h log 

r a — r a 


Stt^coc 

so that there results the spherically symmetric expression 


- 2 


(19) 


C'p(r,t) 


h. 


1 


dTT^eoC 


3ft 


c(t — ir) , c(t 

-^log-7- 

2r c{t 


ir) + r 
ir) — r 


( 20 ) 


3.1 Small distance corrections to the fields and correlators 

Our aim is to determine a numerically treatable presentation of the stochastic 
helds beyond the dipole (i.e., r —>■ 0) approximation of NLl. We switch from SI 
units to Bohr units. This implies that eA is replaced by /?A, with /3 defined in 
(10) and A being the sum (11) with the argument of the square root replaced by 
STrjVk, involving the normalised variables V = V/{ctq)^ and fc = wtq. In these 
units it holds that Tc = 1, while Cs and Cp read 


Cs{r, t) 
Cp{r, t) 


-:^3ft 


27r {t 
3 


27rr2 


3ft 


1 

— ir)^ — ’ 

t — ir , t — ir + f 

-^—\og- -^-= 

2r t — IT — r 



( 21 ) 


to be taken at r = Zar. Clearly, we need only small r, where these functions have 
leading order behaviours 


Cs{r,t) = 


C'p(r,t) = -—3ft 


+ 


(t — ir)^ {t — it) 
1 3r^ 


{t — irY 5(t — ir)^ 


( 22 ) 
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From (13) and (15) we collect the leading behaviour for the A autocorrelation 


CA(r,t;0,0) 


Iso 1 , 3 rr- 2r^l 

TT {t — iTc)^ hTT {t — ITc)^ 


(23) 


Numerically, it is far too cumbersome to deal with the 3 — d sum over k values 
in the dehnition (11). In NLl it was realised for the limit r —>■ 0 that, since the A 
helds are sums of Gaussian random variables, they themselves are Gaussian and 
may just as well be replaced by different Gaussian sums that produce the same 
autocorrelator. We wish to extend this approach to include the r-dependence to 
quadratic order. Hence we write 


q, s) 



dwlF^(a;)C'4(r, 


(24) 


with W(u!) = e which demands 


Cfj{r,q,ui) = {A^{u:,r)AJ{-uJ,q)) = + qtQj - 2(r^ + g^)d,j] 

TV iUTT 

^3 

+T7;-(4^fcafc Sij - fiQj - fjq^) (25) 

iUTT 

Let us recall the solution for the r = q = 0 case of NLl. First, the fre¬ 
quency integral is discretised. For uniformly spaced frequencies Wn = ndu) with 
n = 1,2,3,---, the leading term may arrive from 


A 

A(0,t) = ^ J W (uJn) ( An sinwnt -h B„ cos Writ ), (26) 

n=l 

having independent, real valued Gaussian random variables A^t and B^i for n = 
1, 2, 3, • • • and i = 1, 2, 3, with zero average and unit variance. 

The new task is to extend this for small r and q. Since the autocorrelator has a 
correction of order r^, we want for A linear and bilinear terms. While the leading 
correction must be linear in r or q, it happens that the term of A bilinear in r 
can be determined explicitly. This leads to the component of A 


Ai (Wtj , r) — Ajii ^nik^k T ['^i^k'^k ‘2ciiT ], 


(27) 


for i = 1,2,3 and with sum over k = 1, 2, 3 implied. 
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To consider the terms linear in r, we introdnce 8 real valned matrices, which 
coincide up to a factor with the generators of SU(3), 


/O 1 0 \ 

= ^/3 1 0 0 , 

Vo 0 0 / 

/I 0 0 \ 

A^ = ^/3 0 -1 0 , 

Vo 0 0 / 

/O 0 -1 \ 

A'^ = \/5 0 0 0 , 

Vl 0 0 / 

/O 0 0 \ 

a'^ = \/5 0 0 -1 , 

Vo 1 0 / 


/O -1 
A^ = x/5 1 0 

Vo 0 

/O 0 
A^ = \/3 0 0 

Vl 0 

/O 0 
A® = \/3 0 0 

Vo 1 




(28) 


At each frequency Wn > 0 we express the random bnik in terms of a scale Cn and 
8 independent, complex valued Gaussian random variables Pna, 


8 



a=l 


The P’s are chosen to have averages zero and unit variances in the manner 


Pna = - ipi% {Pi^^Pl^l)={Pi?Jp‘;^l)=SabSnm, (30) 

while by assumption {PnJ) = (PnJ) = {PnJP^^l) = 0. This implies that 

iPnaPmb) = ^SabSnm {u, b = I, ■ ■ ■ , 8), (31) 

while a bit of algebra shows that 


8 

yy Kk^ji = ‘2.{4:6ij6ki — SikSji — SiiSjk ). 

a=l 


Hence the choice 


Cn 


da; u!n 
207r 


reproduces at each n the desired cross terms of r and q. 


(32) 


(33) 


{bnikbnji)rkqi = i^rkQkSij - ngj - rjQi). 


with sums over k and I implied. We end up with the field up to order 


(34) 
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(^n) { (-^?^^ COS COrii H~ -^ni Sin LUyit^ 

n>0 


+ cosujnt + sinu;rii)A?fcrfc (35) 

4\/5 
^2 

ni‘^ COSUJnt + {AnkThTi - 2Anir^) SinWnt]} 


summed over k = 1,2,3 and a = I,-- - ,8. Notice that the Coulomb gauge con¬ 
dition diAi = 0 is satished and that the translation invariant autocorrelator (23) 
comes out right up to order included. 

We now acknowledge that the right hand side of (35) should be taken at f = 
Zar. The electric field Ei{r,t) = —dtAi{r,t) thus follows as 


/dZj 

E^{r,t) = Y, y - ^^^—^W(uJn){{BniSinUJnt - Ani COSUJnt) 

+ sinwnt - cosa;„t)A“fcrfc (36) 

4v5 

-h Z'^a‘^^[{Bnkrkri - 2Bmr^) sinuj^t - (Ankrkn - 2Anir‘^) cosw^t]}. 

As is well known [1,11], the B field is of order Za. Indeed, the field tensor Fij — 
diAj — djAi reads 


Fij = -ZaY W{uj-a) Y + ( 3 ’^) 

n>0 ^ a=2,5,7 

z'^a'^ Y W{LOn)[{Bnirj - Bnjri) COSUnt -h {AniTj - AnjVi) SinWnt], 

n>0 

containing only the antisymmetric A“. The magnetic field components follow as 
Bk = ^SijkFij and the Lorentz force reads (r x B)i = Fijr-j. 

With the electric and magnetic fields expressed as single sums over the fre¬ 
quency, we can now implement this structure in the numerical code. In our OpenCL 
method not much extra time is demanded for handling at each cOn the 16 additional 
random variables /3na'^\ 


4 Numerics 

In our code we numerically solved the ODE in (9) for a charge of Z = 3 using 
a fourth order Runge-Kutta scheme. For our most important run, we cut off the 
frequency in equations (36) and (37) using a ‘moving’ cut-off at 2.5 times the 
Keplerian frequency of the electron. The used code is an enhanced version of the 
code described in [12]. This code was/is highly optimized using OpenCL and it 
is run on a state-of-the art GPU (AMD R9 290X). The present performance is 
equivalent to a CPU-only version of the code running on a 300 core CPU cluster. 
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Here we give a pointwise summary of the main features/caveats and the differences 
to the version described in ref. [12]: 

1) The considered frequencies are Wn = n/N. Due to the presence of a vast 
amount of random variables we had to set N = 10^ instead of = 10® as was 
done in our previous article [12]. Up to t = N the autocorrelator of equation (35) 
reduces to the autocorrelator (23) of the exact A field, defined in equation (11). 
For larger t, we may not simulate a genuine 3D problem anymore. This could be a 
problem since it is shorter than the damping time l//3^ and possibly also shorter 
than other timescales relavant to the electron’s dynamics. 

2) We used the ’classical’ Runge Kutta fourth order integration scheme with 
full energy and angular momentum conservation up to eccentricities of 0.99 and 
higher. We did also test other ODE schemes like several of the Adam-Bashforth 
type and the simpler Euler method. On top of that we did vary the number of 
iterations per orbit from run to run. Thus we can conclude that our results are 
ODE scheme independent, with the biggest numerical error induced due to the 
other points mentioned in this paragraph. 

3) We used 4000 iterations per orbit. The EM-field, though, was only updated 
25-30 times per period of the highest mode with the other points determined by 
higher order interpolation using Lagrange Polynomials. We calculated this EM- 
field 1.5 to 2 times as often with respect to our previous article [12], since higher 
numerical precision was required so to assure that the effects of the relativistic 
corrections, the magnetic field, the positional dependence of the EM-field and the 
spin-orbit coupling were fully accounted for. 

4) As was done in our previous article, we updated the moving cut-off frequency 
when the orbital period of the electron changed more than 20%. This was done to 
minimise the presence of discontinuities in the equation of motion (9). Note that 
this field switching could still induce errors in our solution. 

5) The code spends a very long time to simulate the electron at low energies. 
To speed the code up, we gave the electron a push to increase its energy when it 
became smaller than £ = —1.6. This is justified by the resulting energy distribution 
and the earlier published conjecture [13]. For the present article we chose a scheme 
in which we randomly gave the electron a push either parallel or perpendicular 
to its velocity axis. The advantage of this scheme is that it remains simple and 
induces less bias in the angular momentum at low energies. 

6) We made use of the Runge Kutta RK4 algorithm. This is not the main 
source for our errors, because others issues like the switching of the EM field every 
few hundred steps, the finite cutoff frequency and other points mentioned above, 
produce (much) bigger errors. We did do extensive convergence testing, though, 
even runs with full double precision and 4 times as many steps as used in the 
reported results (but with smaller N values). No statistical difference was detected 
in any alternative simulation. On top of that, we experimented with several Adam 
Bashforth like methods, again no difference was detected at all. We could even have 
used the standard Euler method. Overall, we have more than 30 runs in different 
configurations and with slightly varying parameters, all giving approximately the 
same results. Since we only investigate the statistical behaviour of SED, but do 
not claim to solve it with high precision at a very late times, we are convinced of 
the fact that the reported results give a fair presentation of the situation. 
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5 Results 

We ran the code several times using different values for the cut-off frequency. The 
longer runs took 50 to 200 GPU hours to complete. 

Here we present the results of our most promising run up to the point that 
the electron energy went above £^=-0.05, what we dehne as ionisation. The moving 
cutoff was set at 2.5 times the electron’s Keplerian frequency. In hgures la,lb we 
show the time series for the energy and radius while in figures 2a, 2b we show the 
corresponding probability distributions compared to the Is ground state and an 
earlier published conjecture [13]. This conjecture for the classical distribution of 
the energy E and the angular momentum L reproduces the quantum Is ground 
state results. In table 1 we give the run time of this simulation before ionisation 
and we compare it to the damping timescale of 1//?^. 


property 

value 

duration (s) 

tfotal 

2.05 10'to 

5.510“"" s 

idamp 

4.28 10'’to 

1.15 10“"'" s 

^orbit 

3.2610'’ 


damp 

48 



Table 1: Time and number of orbits for our simulation with 2.5 harmonics and Z = 


lei 






5.0x10* l.OxlO’ 1.5x10’ 2.0x10’ ^ 


5.0x10* 1.0x10’ 1.5x10’ 2.0x10’ ^ 

Fig. 1 a): Time series for the energy, in Bohr units, b): Time series for the radius. 



Fig. 2 a): Histogram for the energy data of Fig. la. The curve is the conjecture P{S^ = 
(4/3|£^|®) exp(—2/|£^|) from refs. [13,12]. b): Histogram for the radius data of Fig. lb. The 
curve is the quantum prediction P = 4r^ exp(—2r). 



















































































H-atom in Stochastic Electrodynamics 


11 


These results suggested that the electron remained stable and had a probability 
distribution closer to the conjecture [12, 13] for a longer time than was presented 
in our previous article [12]. However, we did find out that the time up until the 
electron’s ionisation varied quite drastically between t = 10®fo to t = lO^to for 
different runs using either the old code or the new code. Thus we can conclude that 
there is no statistically significant difference between this run and the idealised run 
in our previous article [12]. Runs with the moving cutoff set at 4.5 or 6.5 times 
the orbital frequency showed ionisation on even shorter timescales. Multiple runs 
with a fixed cutoff set at 1.5 times the Kepler frequency for £ = —1.6 did even 
show ionisation within t = lO^to, though the accuracy of the interpolation scheme 
may be brought into question for the fixed cut-off case, especially when we want 
to account for the higher order effects discussed in this article. This is because for 
high electron energies very strong EM-modes are included, the frequency of which 
is 30 times the Keplerian frequency of the electron, while we expect only the first 
few ‘harmonics’ to have a significant contribution [11, 12]. Therefore the numerical 
error in the value of the EM-field can become as large as 20 percent of the first 
harmonic. 


6 Discussion 

We have extended our previous dynamical simulation for the hydrogen ground 
state in Stochastic Electrodynamics (SED) to include relativistic corrections: the 

kinetic correction and the spin-orbit coupling. Formally these terms are larger 
than the Lamb shift, which should be responsible for effective quantum behaviour 
of SED. While in our previous approach [12] a self ionisation was observed, the 
question is how the relativistic corrections - stronger than the SED corrections 
- impact on this. Hereto the stochastic fields, sums over 3d fc-space, are first 
expressed as frequency sums of stochastic variables that reproduce the same cor¬ 
relations up to the needed order r^. Next, the routine is implemented in OpenCL. 

The answer to the stability question is negative: the relativistic corrections 
have hardly any impact: self ionisation remains present; the relativistic corrections 
just act as small terms to the non-relativistic Hamiltonian. This is at odds with 
a prediction put forward by Boyer [14,15], who argues in favor of a relativistic 
treatment to arrive at the quantum result for the H atom. This is now shown not 
to work out. The continued finding of self-ionisation underlines the correctness of 
leaving out the relativistic corrections at the most basic level of the theory, as is 
indeed routinely done in the field. The observed discrepancy between SED and 
quantum mechanics may also may carry ingredients of the SED description of the 
free particle [12,16]. 

Whether the thus observed self ionisation is the final nail in the coffin of SED, 
or, perhaps, due to approximations made, remains an open question. The fact 
that ionisation does not happen (in the employed window approximation) at early 
times but at some later time is likely due to the fact that the electron has to 
find the niche in which it can ionise and stay there sufficiently long. An issue of 
special interest is the point charge approximation for the electron. We stress in this 
respect that the point particle theory is incapable to deal with the Darwin term 
of Eq. (2), a relativistic correction proportional to a delta-function. In this sense, 
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the employed theory could not have been correct from the start. Whether these 
deficiencies can be simultaneously overcome remains a challenge for the future. 
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